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3.  SUMMARY 

This  report  summarizes  the  results  of  the  study  undertaken  by  the  authors  under  USAF  Grant  n°  FA9550-16- 
1-0081. 

The  work  performed  lies  within  the  context  of  a  wider  collaboration  agreement  undertaken  jointly  by 
UniBo  team  with  a  private  company  and  a  research  institute,  towards  the  development  of  a  MEMS 
prototype  of  a  hydrogen  peroxide  microthruster  delivering  thrust  in  the  mN  range,  and  targeting  the  class 
of  micro-nano  spacecraft. 

Most  of  the  work  consisted  in  theoretical  and  numerical  analysis  of  the  flow  within  the  dissociation 
chamber.  This  is  deemed  to  be  an  important  tool  to  guide  the  design  phase  of  the  microthrusters. 

The  report  is  structured  as  follow: 

•  Introduction  and  motivation  for  the  investigation. 

•  Development  of  the  numerical  model,  including  detailed  discussion  of  the  mathematical  framework 
and  analysis  of  the  results. 

•  Design  activities  in  preparation  of  the  manufacturing  of  the  next  thrusters  prototypes. 

•  Conlcusions. 
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6.  INTRODUCTION 

Many  Space  Agencies  are  currently  pursuing  the  development  of  light,  small  S/C  platforms  for  a  wide  range 
of  missions  (such  as  Earth  observation  and  formation  flying),  and  a  major  critical  point  is  indeed 
represented  by  the  on-board  propulsion  technologies  to  meet  the  orbit  correction/control  and  de-orbiting 
needs.  The  typical  AV  values  needed  by  such  missions  are  incompatible  with  simple  cold-gas  systems  but,  at 
the  same  time,  the  toxicity  hazards  placed  by  hydrazine  or  the  often-prohibitive  power  needs  of  electrical 
propulsion  calls  for  the  search  of  alternative  propellant  solutions.  A  very  promising  technology  is 
represented  by  hydrogen  peroxide  based  monopropellant  systems,  which  offers  the  highest  volumetric 
specific  impulse,  as  compared  to  many  other  propellants. 

Propulsion  systems  in  which  the  release  of  thermal  energy  is  obtained  through  a  catalytic  decomposition, 
need  a  propellant  only  and,  for  this  reason,  are  usually  called  "monopropellant  thrusters". 

Due  to  its  reduced  toxicity,  the  nature  of  the  decomposition  products  and  environmental  friendless, 
hydrogen  peroxide  has  gained  increased  attention.  In  micro/nano  spacecraft  applications  where  cost 
reductions  are  being  emphasized,  eased  procurement  and  handling  procedures  that  such  a  propellant  yield, 
feature  prominently.  On  top  of  that,  a  great  asset  in  miniaturization  may  come  from  a  MEMS-based 
thruster  manufacturing  process,  which  has  been  already  exploited  in  the  past  years  for  cold-gas  micro¬ 
propulsion. 

6.1  Motivations 

One  of  the  challenges  with  microscale  thrusters  is  that,  as  the  size  is  reduced,  the  ratio  of  the  surface  area 
to  volume  increases.  Therefore,  the  heat  losses  (surface  phenomenon)  become  significant  compared  to  the 
heat  generated  by  the  exothermic  reaction  within  the  device  (voume  phenomenon),  and  the  thermal 
energy  available  inside  the  combustion  chamber  may  not  be  sufficient  for  sustaining  the  reaction  and 
vaporization  of  the  liquid  products  and/or  non-reacted  peroxide.  The  effects  of  boundary  layers  also 
become  more  pronounced  at  small  scales  and  can  lead  to  significant  deviations  from  inviscid  behavior. 

For  similar  reasons,  MEMS-based  devices  are  well  suited  to  host  heterogeneous  catalytic  processes,  which 
are  surface  reactions,  due  to  their  intrinsically  high  surface  to  volume  ratio.  However,  the  lack  of  sound 
understanding  of  the  coupling  between  heterogeneous  catalysis,  multiphase  fluid  dynamics,  and  thermal 
phenomena  makes  it  a  difficult  task  to  estimate  the  critical  length  needed  to  obtain  a  complete  propellant 
decomposition,  based  on  simple  considerations  only. 

An  empirical  approach  would  be  the  one  which  scales  the  literature  data  available  for  larger  HTP  thrusters 
down  to  the  MEMS  size.  This  exercise  results  in  an  estimated  length  of  =  2  mm.  At  this  small  scales  the 
Reynolds  number  is  very  low  leading  to  laminar  flow  regardless  of  the  geometry  and  therefore  very  little 
mixing  occurs.  In  other  words,  diffusion  and  conduction  are  the  dominant  modes  of  species  and  heat 
transport.  This  in  turn  suggests  that  direct  application  of  downscaling  is  questionable,  since  in  large  scale 
devices  turbulent  mixing  occur.  It  is  therefore  reasonable  to  expect  the  former  exercise  to  lead  to  a 
considerable  underestimation  of  the  critical  length. 

The  critical  channel  length  must  be  selected  such  that,  for  a  given  mass  flow  rate,  complete  decomposition 
occurs.  If  the  chamber  is  too  short,  incomplete  decomposition  will  result  in  diminished  efficiency.  If  the 
channel  is  too  long  condensation  may  occur. 
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Given  the  discussion  above,  the  design  of  the  thrusters  object  of  this  study  cannot  be  based  solely  on  semi- 
empirical  arguments  rather,  a  comprehensive  numerical  analysis  is  required. 

6.2  Background 

The  background  activity  undertaken  by  the  team  as  a  joint  collaboration  agreement,  lead  to  a  first  design  of 
microthruster  a  few  years  ago:  this  design  forms  the  baseline  for  the  next  generation  also.  The  first 
prototypes  of  H202  micro-thrusters,  manufactured  at  an  external  research  institute  premises,  were  silicon 
MEMS  devices  featuring  a  decomposition  chamber  configured  as  a  dense  pattern  of  staggered  pillars.  The 
great  number  of  pillars  was  aimed  at  maximizing  the  catalyst  surface  per  unit  volume  to  enhance  the 
decomposition.  The  overall  dimensions  of  those  prototypes  were  approximately  18x4xlmm3. 


Figure  6-1:  Staggered  pillar  concept  for  the  decomposition  chamber. 

Experimental  testing  of  those  samples  showed  however  an  insufficient  propellant  conversion,  with  a 
bubbly  gas-liquid  mixture  found  at  the  outlet.  These  results  called  for  a  more  in  depth  investigation  of 
the  underlying  physics,  to  provide  some  guidelines  for  the  design. 

For  the  numerical  analysis  performed  during  this  study,  the  staggered-pillars  layout  was  retained:  indeed, 
as  a  result  of  the  design  optimization  we  expect  to  tailor  the  size,  spacing,  shape  of  pillars,  without 
changing  dramatically  the  overall  concept. 

7.  NUMERICAL  MODEL 
7.1  Literature  Review 

The  design  of  a  reliable  high  test  peroxide  based  monopropellant  systems  is  strictly  related  to  our  capacity 
of  developing  a  numerical  model  capable  of  simulating  the  heterogeneous  catalytic  decomposition  process 
of  liquid  hydrogen  peroxide. 

The  exothermic  reaction  at  the  base  of  the  operation  of  the  device  evolves  as  a  sequence  of  the  three 
following  stages: 

1.  A  first  stage  where  both  peroxide  and  water  exist  in  the  liquid  phase; 

2.  A  second  stage  where  the  aqueous  peroxide  mixture  is  boiling; 

3.  A  third  and  final  stage  where  both  are  in  the  gas  phase. 

In  all  three  stages,  peroxide  and  water  coexist  with  gaseous  oxygen. 

It  is  evident  as  a  full  simulation  of  the  decomposition  reaction  is  extremely  challenging  and  complicated  by 
the  multiphase  nature  of  the  flow  and  by  the  coupling  of  heat  and  mass  transfer. 
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There  is  a  quite  extensive  literature  on  numerical  modelling  for  the  decomposition  of  hydrogen  peroxide  in 
thrusters  and/or  micro-reactors  in  general.  As  a  matter  of  fact,  this  topic  cannot  yet  be  considered  fully 
covered,  since,  whether  to  a  smaller  or  to  a  greater  extent,  some  simplifying  assumptions  shall  be  made  if 
one  wants  to  develop  a  tractable/solvable  numerical  model. 

A  common  simplified  approach  is  the  one  which  models  one  single  phase  at  a  time.  For  example,  in  [1]  the 
decomposition  of  the  vapor  phase  was  considered  in  a  2D  model  of  a  microreactor,  with  the  chemical 
reactions  modelled  using  a  simplified  formulation  with  one  species  concentration  equation  (typically  used 
for  dilute  solutions),  and  no  interaction  with  the  solid  domain  considered. 

Alternatively,  it  is  possible  to  consider  what  happens  in  a  pre-vaporization  condition  to  evaluate  the  effect 
of  the  monopropellant  in  its  liquid  form  and  the  heat  generation  associated  with  the  decomposition,  as 
done  in  [2] ,  where  the  interaction  with  the  substrate  was  still  absent. 

In  [3]  ,  a  3D  simulation  for  the  decomposition  of  a  50%  aqueous  peroxide  solution  within  a  packed  bed 
microreactor  is  presented,  where  evaporation  is  partially  accounted  for,  but  only  for  the  water  species. 
Again,  the  species  conservation  was  implemented  with  a  simplified  single  equation  for  molar 
concentration. 

Numerical  models  which  encompasses  all  the  phases  have  been  developed  in  the  past,  but  are  limited  to 
ID  studies,  as  in  [4] ,  [5]  . 

As  a  natural  complement  to  the  work  performed  in  the  cited  references,  in  the  present  study  we  sought  for 
a  unique  model  capable  of  modelling  the  decomposition  process  both  in  the  pre-  and  post-vaporization 
phases  of  the  liquid  reactants  and  products. 

7.2  Development  of  an  evaporation  model  for  the  decomposition  of  H202  suitable  for  FEM 
analysis 

The  literature  review  showed  that  three  main  aspects  are  overlooked  in  the  models  developed  so  far: 

1-  A  sound  treatment  of  the  chemical  species  conservation,  usually  relegated  to  a  single  species 
equation,  more  suited  for  diluted  solutions. 

2-  A  theoretical  model  which  jointly  accounts  for  the  dissociation  of  the  propellant  and  its  evaporation 
due  to  the  temperature  increase,  to  be  integrated  in  the  overall  FEM  model. 

3-  Effect  of  the  conjugate  thermal  problem  of  the  fluid  stream  and  the  silicon  substrate. 

This  last  point  is  of  great  importance  due  to  the  high  thermal  conductivity  of  the  silicon  substrate,  which 
enhances  the  transport  of  the  generated  heat  of  reaction  all  along  the  chamber  length.  It  is  worth  to  be 
noted  that  to  capture  all  the  relevant  phenomena,  especially  point  3,  a  3D  model  is  required. 

To  this  end,  a  commercial  FEM  code  was  employed  (COMSOL  Multiphysics)  which  allows  the  coupling  of 
several  physical  interfaces.  The  global  numerical  model  can  then  be  thought  as  the  sum  of  several  PDE 
solved  within  the  two  domains,  fluid  and  solid,  according  to  the  following: 

a)  FLUID 

1.  Species  transport 

2.  Mixture  continuity 

3.  Momentum  transport 

4.  Energy  transport  (including  phase  change) 

b)  SOLID 

1.  Energy  transport 
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To  limit  the  computational  burden,  the  geometric  domain  is  restricted  to  only  one  elementary  slit  (see 
Figure  7-1)  of  the  combustion  chamber,  without  detailed  modelling  of  the  inlet  and  outlet  regions.  This 
simplification  retains  nevertheless  the  fundamental  structure  of  the  current  design  for  the  combustion 
chamber,  implemented  as  a  dense  pattern  of  staggered  pillars.  Main  geometric  dimensions  are  collected  in 
Table  7-1. 

The  mesh  is  composed  by  slightly  more  than  1  million  of  elements,  prism  in  the  fluid  domain  and 
tetrahedral  in  the  solid  domain.  As  apparent  from  Figure  7-2,  a  finer  mesh  was  set  up  for  the  fluid  domain, 
while  a  much  coarser  discretization  was  chosen  for  the  silicon  substrate. 

Due  to  the  chosen  approach  of  modeling  only  an  elementary  slit,  an  extensive  use  has  been  made  of  the 
"symmetric  boundary"  feature,  which  is  common  to  many  CFD  software  for  exploiting  the  planes  of 
symmetry  of  the  geometry  under  study.  By  referring  to  Figure  7-1,  the  symmetry  planes  are  the  two  lateral 
x-z  sides,  plus  the  top  side,  x-y,  since  only  one-half  of  the  channel  height  is  modelled. 


Description 

Size  [m] 

Length 

9e-3 

Width 

le-4 

Height 

2e-4' 

Substrate  thickness 

4e-4 

Pillar  size  (/x)2 

2.6e-4 

Pillar  size  (/y) 

0.75e-4 

Pillar  spacing  (Ax) 

le-4 

Pillar  spacing  (Ay) 

0.25e-4 

Table  7-1:  Summary  of  the  main  dimensions  for  the  baseline 
geometrical  model. 


Figure  7-1:  Isometric  view  of  a  domain  modelled  with  equispaced,  diamond  pillars 


1  Half-heigth  of  the  fluid  domain  (simmetry  condition) 

2  Length  of  rhombus  diagonals. 
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Figure  7-2:  Detailed  view  of  the  inlet  region:  geometry  (left)  and  mesh  (right) 


7.3  Fluid  domain 

The  mathematical  model  for  the  chemically  reacting  flow  dynamics  is  adapted  from  [6]  ,  where 
conservation  equations  are  formulated  for  the  species  mass  fractions. 

7.3.1  Species  transport 

Transport  of  chemical  species  is  subjected  to  the  following  equation  for  the  mass  fraction  Yk,  describing  the 
conservation  of  mass: 


P 


DJk 

Dt 


-V  ■  jfc  +  rkMMk 


(7-1) 


where  ikis  the  diffusive  mass  flux  and  —  =  rr  +  u  ■  V  stands  for  substantial  derivative. 

,K  Dt  dt 

There  exist  several  theories  for  describing  multi-component  diffusion.  Here  we  have  chosen  a  mixture- 
averaged  formulation,  which  is  not  the  most  rigorous  implementation,  still  it  is  sufficient  for  our  purpose, 
which  is  to  perform  an  "order-of-magnitude"  evaluation  of  the  impact  of  diffusion  coefficients  over  the 
combustion  chamber  efficiency.  In  such  formulation,  the  diffusive  flux  reads  as: 


jk  =  pDr(vYk+ivMM)  (7-2) 

where  MM  and  D™  are,  respectively  the  mixture-averaged  molar  mass  and  mixture-averaged  diffusion 
coefficient  for  species  k3.  Note  that,  to  ensure  total  mass  conservation,  the  equation  above  is  solved  for  2 
species  only,  and  the  mass  fraction  of  the  third  one  is  computed  as  the  complement  to  1.  In  the  simulations 
performed,  the  specie  solved  for  through  mass  conservation  is  H202. 

Note  that  p  is  the  average  mixture  density,  and  it  will  be  better  defined  when  describing  the  energy 
transport  equation  and  the  evaporation  model. 


3  Full  expressions  for  MM  and  D™  can  be  found,  e.g.  in  [6] . 
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Summing  up  the  contribution  of  the  different  species  leads  to  the  standard  continuity  equation  for  single¬ 
phase,  non-reacting  fluid,  since  the  sum  of  all  diffuses  fluxes  as  well  as  all  reaction  rates  must  be  zero  by 
definition: 


(7-3) 


7.3.2  Momentum  transport 

The  momentum  transport  is  described  by  the  well-known  Navier-Stokes  equation: 


—pi  +  /z(Vu  +  Vur)  —  -/z(V  ■  u )/ 


(7-4) 


Where  u  is  the  velocity  field,  p  the  pressure  field  and  p  the  dynamic  viscosity  of  the  mixture.  This 
implementation  relies  on  the  assumption  of  low-Mach  number  flow,  which  is  certainly  satisfied  within  the 
combustion  chamber. 

7.3.3  Energy  transport 

The  equation  for  the  transport  of  enthalpy  (on  a  mass  basis,  i.e.  [J/kg])  in  a  homogeneous  chemically 
reacting  flow  reads  as: 


DH  Dp  V-1 

PDt=Dt  +  V'('CV7’)-Z,V'(Hk'ik)  +  0  (7'5) 

k= 1 

It  highlights  the  various  mechanisms  for  enthalpy  transport,  i.e.  pressure  work,  Dp/Dt,  thermal  conduction 
V  ■  (kVT),  enthalpy  variation  due  to  diffusive  transport  of  species  Sk=i  ^  '  (Mk  '  jk)>  and  viscous  dissipation 
<t>  (whose  explicit  expression  may  be  found  in  standard  textbooks).  In  the  present  application,  (low 
pressure,  low  velocity  flow)  the  viscous  dissipation  and  pressure  work  terms  can  be  safely  neglected. 
Furthermore,  it  is  often  convenient  to  express  the  enthalpy  transport  equation  in  terms  of  temperature 
only,  which  leads  to: 

K  K 

DT  Dp  v  v 

pCp  Dt  =  m  +  v ' (KVr)  “  L  Cp'kjk  ‘ (KVr)  “  Z, HkhMk  +  0  (7'6) 

k= 1  k= 1 

In  general,  for  a  gaseous  mixture,  the  mixture  averaged  properties  are  evaluated  as: 


=  X  Cp'kYk  ’  H  =  X  HkYk  1  P  = 


k= 1 


k= 1 


v  k  Xk 
Lk=1Pk 


(7-7) 
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l.e.  the  specific  heat  and  enthalpy  are  weighted  averages  of  the  respective  species  coefficient  weighted 
through  mass  fractions,  while  the  density  is  a  weighted  average  of  the  volume  fractions.  Later,  these 
expressions  will  be  generalized  to  account  for  the  phase  change. 

The  enthalpy  of  each  species  has  been  assumed  independent  on  the  pressure  (strictly  true  for  perfect 
gases,  approximately  true  for  real  gases  and  liquids  at  low  pressure): 

dHk  cpkdT 

The  temperature  formulation  brings  also  naturally  out  the  concept  of  enthalpy  of  reaction,  through  the 
introduction  of  the  source  term  £jL1/4rfeMMfe  which  depends  on  the  reaction  rate  rk  [mol/m3s]. 

Note  that  the  expression  for  the  mixture  density  is  in  accordance  with  the  usual  volume  fraction 
formulation  for  immiscible  fluids,  as  it  is  the  case  when  dealing  with  different  phases. 

For  the  case  under  study,  where  the  catalytic  dissociation  of  hydrogen  peroxide  occurs  at  the  catalytic 
surface,  the  enthalpy  of  reaction  enters  the  system  as  a  boundary  heat  flux  rather  than  a  volumetric  source 
(and  its  dimensions  change  to  [mol/(m2s)]  accordingly). 

This  formulation,  although  general,  does  not  directly  accommodate  phase  changing  fluids.  However,  the 
vaporization  process  can  be  handled  thanks  to  a  technique  that  is  known  in  the  literature  as  the  "equivalent 
heat  capacity  method",  originally  conceived  for  melting  solids,  and  further  developed  in  this  framework  to 
the  problem  at  hand. 

Essentially  this  method  rather  than  employing  the  common  latent  heat  of  vaporization  Lvap  as  an  additive 
term  in  the  energy  balance  equation  when  a  material  reaches  its  phase  change  temperature  Tpc,  prescribes 
that  the  phase  change  occurs  in  a  temperature  interval  between  Tpc  -  AT/2  and  Tpc  +  AT/2.  In  this  interval, 
the  following  assumptions  are  made: 

1.  A  smooth  function  of  temperature,  6,  is  defined,  which  is  equal  to  1  before  Tpc  -  AT/2  and  to  0  after 
Tpc  +  AT/2.  It  corresponds  to  the  volume  fraction  of  phase  1  within  the  mixture 

2.  The  density,  p,  and  the  specific  enthalpy,  H  [J/kg]  of  the  mixture  are  expressed  as  a  convex 
combination  of  the  respective  values  for  the  phase  1  (before  evaporation  starts)  and  phase  2  (fully 
evaporated  mixture). 

Assumption  2  results  into  the  following  equations: 

p  =  ePl  +  (1  -  0)p2  (7-8) 

pH  =  dp^  +  (1  -  6)p2H2  (7-9) 


The  same  combination  is  also  used  for  computing  the  mixture  viscosity  p,  necessary  for  solving  momentum 
transport: 


p  =  6p1  +  (1  -  0)p2  (7-10) 

Note  that  mass  and  species  convervation  are  ensured  by  using  the  expression  in  Eq.  (7-8)  for  the  density 
used  to  solve  both  Eq.  (7-1)  and  (7-3). 
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Therefore,  the  method  basically  provides  a  constitutive  equation  for  the  variation  of  the  volume  fraction 
within  the  evaporating  mixture  as  a  function  of  temperature  only.  Note  that,  however,  a  dependency  on 
the  chemical  composition  is  still  retained,  since  the  physical  properties  of  the  two  states,  1  and  2,  at  the 
beginning  and  end  of  the  evaporation  interval,  are  indeed  functions  of  the  local  composition  through  Eqs. 
(7-7):  depending  on  the  inlet  propellant  concentration,  and  boundary  conditions,  we  may  expect  that  the 
7pC  -  AT/2  and  Tpc  +  AT/2  temperature  levels  may  be  reached  at  different  degree  of  propellant  dissociation, 
and  this  is  fully  accounted  for  in  the  definition  of  states  1  and  2. 

Differentiating  the  enthalpy  w.r.t.  temperature  to  find  the  specific  heat  capacity,  we  have: 


(7-11) 


Which  is  composed  by  a  first  term,  i.e.  a  volume  averaged  specific  heat  coefficients,  and  a  second  term, 

d0 

proportional  to  the  assumed—,  which  depends  on  the  difference  between  the  enthalpies  of  the  two 

phases,  /-/^  and  H2.  This  last  can  be  interpreted  as  the  enthalpy  of  vaporization  which  is  originated  by  the 
different  enthalpies  of  formation  between  the  liquid  and  vapor  phases  and  is  distributed  over  the  entire 
temperature  range. 

The  specific  enthalpies  for  the  two  phases  are  computed  as: 


Hi  ~  Hq2 Yo2  +  hh2o2(1)Yh2o2  +  HH2o(1)Yh2o 


H2  -  Hq2 yo2  +  hh2o2(v)Yh2o2  +  Hh2o(v)Yh2o 


Where  Ht  =  H°  +  f*0  cp  idT 

The  boiling  point  is  in  principle  a  function  of  the  concentration  of  the  saturated  liquid,  ranging  from  100°C 
(pure  water)  to  150°C  (pure  peroxide),  see  Figure  7-3.  Therefore,  one  may  expect  the  boiling,  within  the 
chemically  reacting  stream  with  variable  composition,  to  be  a  continuous  process  inside  this  interval:  this  is 
the  pivoting  assumption  in  our  method,  which  yields  to  Tpc  =  125°C  and  AT  =  50°C. 
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HI  20  2  Concent  rat  ion,  wt.% 


Figure  7-3:  Plots  of  Boiling  and  freezing  points  of  hydrogen  peroxide  aqueous  solutions  as  a  function  of  weight  fraction 


In  our  case,  we  should  better  talk  of  pseudo-phases  since,  actually,  what  we  can  identify  as  "phase  1"  is 
rather  a  combination  of  liquid  H202-H20  mixture  and  gaseous  02.  Phase  2  is  a  fully  gaseous  mixture  (02  plus 
H202-H20  vapors)  which  is  treated  as  a  perfect  gas. 

The  density  of  phase-1  is  computed  in  the  assumption  of  immiscible  fluids  sharing  the  same  velocity  field 
as: 


1 

Pl  ~~  +  (7-12) 

Pi  Po2 

Where  pt  is  the  density  of  the  liquid  H202-H20  mixture,  known  through  fitting  functions  of  the  temperature 

MMq0V 

and  wt%  available  from  literature,  and  the  oxygen  density  is  computed  as:  p0  — - — . 

^  RT 

MMap 

As  concerning  density  of  phase  2,  from  ideal  gas  law  pg  —  ,  with  the  mixture  molar  mass  MMg  defined 

as: 


MMmir  = 


+  • 


\MMy,  '  MM, 


+ 


MM0 


(7-13) 


7.3.4  Reaction  kinetics 

Reaction  kinetics  has  been  assumed  first  order  in  the  molar  concentration  CH2o2  [mol/m3]: 

1 

i~H2o2  ~  —k'' '  Ch2o2  i~H2o  —  ~^h2o2  fo2  —  —  2^h2°2  [mol/m  s]  (7-14) 

It  is  common  to  assume  for  the  kinetic  constant,  k",  a  temperature  dependency  through  the  so-called 
Arrhenius  exponential  form,  parametrized  with  the  activation  energy  Ea  [J/mol]  and  pre-exponential  factor 
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As  [m/s].  The  former  is  a  property  of  the  catalytic  material;  the  latter  is  more  dependent  on  the  specific 
integration  and  deposition  process  of  the  catalyst  over  the  support. 

Note  that  k"  is  perhaps  the  most  important  parameter  dictating  the  effectiveness  of  the  catalytic 
decomposition  within  the  dissociation  chamber  (the  highest  the  better).  In  the  numerical  analysis,  it  will  be 
assumed  as  an  independent  parameter  in  the  simulations:  clearly  one  must  choose  reasonable  values 
within  a  certain  range.  However,  it  is  not  possible  to  know  a-priori  the  "true"  value  k" ,  that  one  will 
encounter  during  experimental  activities,  since  it  depends  on  many  factors  such  as  the  nature  of  the 
chosen  catalyst  material,  but  also  the  technological  process  used  for  the  deposition  of  the  catalyst  on  the 
pillars'  surface. 

The  kinetic  experiments  performed  at  our  Lab  on  two  commercial  granular  catalysts  showed  that  in  the 
range  of  temperature  23  -  50°C  the  activation  energy  is  in  the  order  of  50  kJ/mol,  which  is  consistent  with 
available  literature  data.  A  direct  application  of  the  kinetic  constant  experimentally  determined  is 
questionable,  since,  as  said  above,  it  is  dependent  on  the  way  the  catalytic  material  is  integrated, 
nevertheless  this  is  a  good  reference  point. 

Much  of  literature  data  reports  the  kinetic  constant  on  a  volumetric  [1/s]  basis  (more  useful  for 
homogeneous  catalysis),  rather  than  a  surface  [m/s]  basis,  as  it  is  required  by  the  FEM  code 
implementation.  The  few  available  reporting  pre-exponential  factor  in  m/s  are  quite  spread  in  the  range  of 
103  to  10s  m/s,  yielding  to  k  in  the  range  [10  4  - 10'1]  m/s  for  temperatures  up  to  =  430  K;  these  values  are 
expected  to  increase  further  if  vaporization  completes  and  fluid  temperature  grows  beyond  such  value.  An 
extensive  compilation  of  literature  data  can  be  found  in  [5]  and  references  therein. 

For  the  simulations  presented  here,  values  of  k  were  drawn  from  the  above  range,  however,  a  simplifying 
assumption  has  been  made.  Since  the  Arrhenius  exponential  dependency  was  found  to  prevent 
convergence  of  the  numerical  solver,  we  chose  to  model  the  reaction  kinetics  assuming  independency  over 
the  temperature.  This  way,  when  drawing  considerations  starting  from  the  results  of  the  simulations  we 
should  regard  the  prescribed  constant  k  as  an  "average"  value  of  the  kinetic  constant  to  be  achieved  under 
experimental  conditions  across  the  reactor  length  in  the  estimated  temperature  range. 

Note  that  this  strong  assumption  it  is  actually  justified  for  the  present  study:  as  will  be  apparent  when 
analyzing  the  simulation  results,  the  presence  of  the  highly  conductive  silicon  substrate  limits  the  fluid 
temperature  jump  between  the  beginning  end  the  end  of  the  catalytic  pattern  to  few  tens  of  degrees. 
Therefore,  only  relatively  small  variations  of  the  kinetic  constant  are  expected  due  to  temperature 
variation. 

7.4  Solid  domain 

The  energy  equation  in  the  solid  domain  (silicon  pillars  and  substrate)  reads  as: 

dT 

PCp  —  =  V-(kVT)  +  Q  (7-15) 

Where  Q  [W/m3]  is  any  heat  source  placed  inside  the  material,  such  as  due  to  joule  heating,  which  is 
applicable  in  case  of  a  resistive  element  placed  in  the  substrate.  The  fluid  and  solid  domain  are  coupled  by 
means  of  energy  exchange  through  the  common  boundaries  by  prescribing  the  same  temperature  at  the 
interface  (i.e.  no  solid-fluid  temperature  slip):  this  feature  is  handled  automatically  by  the  FEM  software. 
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The  exterior  solid  boundaries  are  assumed  to  be  coupled  with  the  surrounding  ambient  through  natural 
convection,  by  letting: 


q  =  Hnc(Ta-T)[W/m 2]  (7-16) 

With  T0  =  298K  and  Hnc  =  40  W/m2  is  a  value  representative  of  natural  convection  through  air. 

7.5  Set-up  for  the  simulations 

Simulations  have  been  performed  with  an  inlet  mass  flux  of  0.1  mg/s.  Outlet  pressure  is  at  atmospheric 
condition.  Note  that,  since  only  one  elementary  pattern  of  the  staggered  pillars  has  been  simulated,  to  get 
an  indication  of  the  actual  mass  flow  through  the  whole  device  one  should  multiply  by  the  total  number  of 
pillars  rows  and  accounting,  at  the  same  time,  for  the  momentum  losses  due  to  the  chamber  lateral  walls. 
Overall,  one  should  expect  the  design  mass  flow  rate  of  the  combustion  chamber  being  a  fraction  of  the 
CFD  single-row  mass  flow  multiplied  by  the  number  of  pillars  rows. 

Several  attempts  were  made  to  enhance  convergence:  as  a  matter  of  fact,  the  formidably  coupled  system 
of  equations  which  create  the  whole  numerical  model  pose  quite  a  challenge  to  the  converge  of  the 
solution. 

The  parameters  used  for  the  baseline  test  case  are  summarized  in  Table  7-3. 

Three  modes  for  running  the  simulations  were  considered: 

1-  Steady  state  simulations  using  a  continuation  technique  for  the  k  parameter.  This  consists  of 
multiple  steady  state  simulations  run  at  increasingly  higher  values  of  k  (called  a  continuation 
parameter)  spanning  from  extremely  low  values  (10'7)  up  to  1-10'3  m/s.  The  solution  for  each  value 
of  k  was  then  used  as  initial  guess  for  the  next  iteration  of  the  continuation  loop. 

2-  Transient  simulations  with  a  constant  value  of  k  and  initial  conditions  taken  from  a  previously 
converged  solutions  at  lower  k  .  These  tests  were  done  basically  in  an  attempt  to  enhance 
convergence  of  the  solutions. 

3-  Transient  simulations;  a  smooth  profile  was  assigned  to  the  k  parameter  up  to  its  regime  value. 
These  tests  give  an  indication  of  a  start-up  phase  of  the  device. 

In  general,  steady  state  simulations  were  achieved  only  up  to  a  value  of  k"  =  2e-3;  to  simulate  higher  kinetic 
rates,  it  was  then  mandatory  to  switch  to  transient  simulations,  whether  using  approach  2  or  3.  These, 
however,  increased  the  computational  burden  even  more. 

Despite  many  "tuning"  efforts  to  enhance  convergence  of  the  simulations,  to  date  reaching  complete 
evaporation  within  one  simulation  was  not  possible,  due  to  errors  encountered  by  the  solver.  Nevertheless, 
a  very  good  insight  could  be  achieved  by  the  simulations,  which  allowed  to  simulate  the  evaporation  within 
the  dissociation  chamber  to  a  degree  higher  than  90%. 
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7.5.1  Summary  of  Boundary  Conditions  and  main  parameters  setting 

In  the  following  table,  a  summary  of  the  BC's  applied  to  the  domain  for  the  different  physics  are  collected 


Boundary 

Boundary  Condition 

Equation 

Mass  Balance 

Inlet 

Mass  fraction 

Y,  -  Y* 

Outlet 

No  diffusive  flux 

-n  •  pD(mVY;  =  0 

Reactor  walls 

Prescribed  total  flux 

-n  ■  p(D,mVY;  +  uYt)  =  f,  Af, 

Bottom  walls 

Null  total  flux 

-n  ■  p(D,mVY;  +  uYt)  =  0 

Symmetry  walls 

Momentum  Balance 

Null  total  flux 

-n  ■  p(D,mVY;  +  uYt)  =  0 

Inlet 

Mass  flow  rate 

—  1  p  (n  ■  u)dS  =  m 

Jd  n 

Outlet 

Pressure 

2 

-pi  +  p(Vu  +  Vur)  -  - p(V  ■  u) 

=  ~Po 

Reactor  Walls 

No  slip 

u  =  0 

Bottom  Walls 

No  slip 

u  =  0 

Symmetry  walls 

Energy  Balance,  Fluid  domain 

Symmetry 

nu  =  0 

Inlet 

Interface  with  isothermal 
domain  at  T=T0 

n  ■  ZcVT  =  pux[H(T0)  -  H(T )] 

Oudet 

No  conductive  flux 

n ■ kV T  = 0 

Reactor  walls 

Boundary  Heat  Source 

n  ■  /cVT  =  Qreac 

Symmetry  walls 

Energy  Balance,  Solid  domain 

No  conductive  flux 

n ■ fcVT  = 0 

Reactor  walls 

Temperature  continuity  with 
fluid 

T-Tfi 

External  walls 

Convective  flux 

n  ■  /cVT  —  Hnat  (Text  ~  T) 

Symmetry  walls 

No  conductive  flux 

n  ■  /cVT  =  0 

Table  7-2:  Summary  of  boundary  conditions  for  the  CFD  decomposition  chamber  analysis 


Symbol 

Value  [units] 

Description 

k" 

from  le-7  to  5e-3  [m/s] 

Average  kinetic  constant 

in 

le-7  [kg/s] 

Inlet  mass  flow  rate 

D? 

le-6  [m/s] 

Mixture  average  diffusion  coefficient  (common  to  all 
species) 

sp 

le-4  [m] 

Pillars  spacing  (channel  width) 

Hnat 

40  [W/ (m2K)] 

Coefficient  of  natural  convection  heat  transfer 

PJ 

0  [W] 

Heat  supply  by  joule  effect 

Tpc 

398.15  [K] 

Assumed  phase  change  temperature 

AT 

50  pq 

Assumed  width  of  phase  change  T  interval 

Table  7-3:  Summary  of  baseline  test  cased  parameters  considered  for  the  CFD  decomposition  chamber  analysis 
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7.6  Results 

7.6.1  Numerical  consistency  of  the  solution 

It  is  known  that  Finite  Element  methods  do  not  exactly  solve  the  conservation  equations  within  the  domain, 
but  only  approximately.  Therefore,  it  is  useful  to  address  to  which  extent,  on  a  global  domain  scale,  the 
conservation  equations  are  satisfied.  To  this  end,  the  following  two  tests  have  been  set  up  which  focus  on 
the  mass  and  energy  conservation: 


1 


puxY  H2o2dAin 


-II 


puxYHj0?dA 


out 


-si : 


i'H2o2MH2o2dA 


cat 


Pout  ~  Pin  *  ||  P^-x^(Tout)dAOUf  +  ||  p  itx  [// (T )  H(T0)]dAi 

JJa,„ 


Aout 

+ 


II  Hnat(Text  T)dAext  —  II  pUxH(Tin)dAin  +  ||  QreacdAcat 
JJ  a _ *  JJa:~  JJa-„+ 


(7-17) 


(7-18) 


which  are  applicable  only  for  simulations  at  steady  state. 

The  first  equation  above  dictates  that  the  difference  of  the  outgoing  peroxide  mass  flux  to  the  ingoing  flux 
must  equate  the  peroxide  mass  flux  reacting  at  the  catalytic  surfaces  (Acat).  This  condition  also 
automatically  guarantees  that  the  enthalpy  of  reaction  is  consistent  with  the  degree  of  conversion 
predicted  by  the  simulation.  The  second  equation  requires  that  the  energy  entering  the  computational 
domain  (i.e.  enthalpy  of  the  inlet  stream  plus  heat  of  reaction)  equates  the  energy  leaving  the  system 
through  three  mechanisms.  These  are:  1)  the  enthalpy  of  the  outlet  stream,  2)  the  losses  through  external 
silicon  surfaces  towards  the  ambient,  and  3)  as  the  enthalpy  lost  to  heat  up  the  propellant  from  its 
isothermal  reservoir  condition  at  Ta  up  to  the  computed  Tin  at  the  inlet. 

In  all  steady-state  simulations,  the  above  equality  was  satisfied  with  an  accuracy  within  ~1%. 


7.6.2  Flow  topology 

Flow  results  are  presented  here  through  a  series  of  3D  plots  of  the  relevant  flow  properties:  velocity, 
density,  temperature,  hydrogen  peroxide  mass  fraction,  evaporation  fraction.  Plots  are  shown  for  values  of 
k  =  [le-4,  le-3,  4.5e-3]  m/s.  Apart  from  the  first  value  for  which  a  direct  steady  state  solution  was  feasible, 
the  remaining  test  cases  were  in  MODE#2:  the  plots  shown  correspond  to  the  last  time  step. 

In  the  following,  some  considerations  emerging  from  inspection  of  pictures  in  Figure  7-4,  Figure  7-5  and 
Figure  7-6  are  briefly  listed: 

•  Velocity  and  Reynolds  number  are  increasing  considerably  on  the  streamwise  direction:  this  is  a 
direct  consequence  of  the  decreasing  mixture  density  spanning  three  orders  of  magnitude  from  * 
103  down  to  10°  kg/m3. 

•  Pressure  drop  across  the  pillars  modelled  chamber  length  is  on  the  order  of  3-4  kPa. 

•  Temperature  variations  across  the  pillars  are  quite  limited,  when  compared  to  the  predictions 
made  using  simplified  ID  models  based  on  plug-flow  theory.  This  is  due  to  the  presence  of  a  highly 
conductive  silicon  substrate,  which  makes  any  analysis  neglecting  the  fluid-solid  conjugate  heat 
transfer  problem  hardly  applicable. 
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•  For  the  same  reason,  we  can  notice  that  the  vaporization  is  quite  well  spread  across  the  reactor 
length,  and  not  concentrated  in  a  thin  region. 

•  By  looking  at  the  plots  of  the  peroxide  mass  fraction,  one  notice  that  higher  decomposition  is 
found  close  to  the  reacting  pillars,  which  is  to  be  expected  since  the  reaction  is  not  volumetric  but 
a  surface  phenomenon.  As  a  matter  of  fact,  a  strong  role  is  played  by  mass  diffusivity  ( D ™),  due  to 
the  lack  of  turbulent  flow  mixing  inherent  to  the  low-Re  micro-scale  regime. 


k  =  le-4  m/s 


yh2o2 


Temperature  [K] 


Pressure  [Pa] 


A  2057 
x!0J 


Velocity  magnitude  [m/s] 

A  0.80 


Evaporation  fraction 

A  0.03 
xl  o'1 


Re 


k  =  le-3  m/s 


yh2o2 


Temperature  [K] 
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Pressure  [Pa] 


Evaporation  fraction 
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k  =  5e-3  m/s 


'h2o2 


Temperature  [K] 


Pressure  [Pa] 


Evaporation  fraction 


Velocity  magnitude  [m/s] 


Re 
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Figure  7-7:  Sample  close-up  views  of  the  peroxide  mass  fraction  (left)  and  stramwise  velocity  (right) 


7.6.3  Combustion  chamber  performance 

To  evaluate  the  performance  of  the  combustion  chamber,  we  focus  onto  relevant  fluid  properties  averaged 
at  the  outlet  surface,  such  as  propellant  percentage  conversion4,  temperature  and  vaporization  fraction. 
Results  will  be  reported  for  Mode#3  simulations:  overall,  this  set-up,  was  the  one  allowing  to  get  to  a  result 
with  the  lowest  solution  time.  The  highest  value  of  k"  reaching  convergence  were  between  4.2-10'3  and 
4.5-10’3  m/s.  Attempts  to  get  beyond  those  values  resulted  into  solver  errors. 

.  Different  settling  time  for  the  kinetic  constant  raising  phase  have  been  used  to  explore  possible 
differences  in  the  transient  phases,  when  using  steeper  or  smoother  transients. 

These  quantities  are  collected  for  the  various  simulations  in  the  tables  below.  For  other  relevant  parameter 
values  employed  for  the  simulations,  refer  to  Table  7-2,  Table  7-3. 


Settling  Time  for  k”:  14  s 

Kinetic 

constant 

K  [m/s] 

Outlet 

peroxide  mass 
fraction  [-] 

Propellant 

Temperature 
conversion  r 

[%] 

Evaporation 
fraction  [%] 

le-3 

0.67 

22  405 

73 

4.5e-3 

0.58 

34  412 

92 

Table  7-4:  Results  of  the  simulations  for  different  values  of  the  kinetic  constant,  diamond 

pillars,  uniform  spacing. 


7.6.4  Effect  of  change  in  pillars  base  geometry 

The  effect  of  a  changing  the  pillars'  base  geometry  was  preliminary  addressed,  by  testing  a  square  base,  in 
opposition  to  the  diamond  baseline  configuration.  The  side  of  the  square  base  was  chosen  as  to  maintain 
an  equal  Surface/Volume  (S/V)  ratio,  i.e.  /x=/y=0.140  mm,  according  to  the  nomenclature  of  Table  7-1.  The 
spacing  between  adjacent  pillars  was  also  kept  fixed.  Results  are  summarized  in  Table  7-5 


4  Defined  as  (^2o2_in  -  W„  J/^-m  ‘  100- 
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Settling  Time  for  k”\ 

12  s 

Kinetic 

constant 
fc"  [m/s] 

Outlet 

peroxide  mass 
fraction  [-] 

Propellant 

conversion 

[%] 

Temperature 

[K] 

Evaporation 
fraction  [%] 

le-3 

0.67 

23 

405 

74 

4.2e-3 

0.58 

34 

412 

92 

Table  7-5:  Results  of  the  simulations  for  different  values  of  the  kinetic  constant,  square  pillars, 

uniform  spacing. 


7.6.5  Effect  of  variable  pillars'  spacing 

The  effect  of  variable  spacing  of  the  pillars  was  also  addressed.  In  particular,  a  negative  gradient  along  the 
streamwise  direction  has  been  explored:  the  rationale  for  this  is  that  it  might  be  convenient  to  have  higher 
S/V  ratios  when  the  propellant  concentration  is  low,  and  lower  S/V  when  concentration  is  high,  in  order  to 
have  a  more  uniform  reaction  rate5  along  the  chamber.  The  geometry  simulated  started  from  an  initial 
spacing  which  was  doubled  that  of  the  baseline  geometry  at  the  inlet,  which  decreases  linearly  the  baseline 
value  at  the  outlet. 


Settling  Time  for  k” 

:  6  s 

Kinetic 

constant 
fc"  [m/s] 

Outlet  Propellant 

peroxide  mass  conversion 

fraction  [-]  [%] 

Temperature 

[K] 

Evaporation 
fraction  [%] 

le-3 

0.56 

36 

364 

0 

4.5e-3 

0.59 

33 

411 

91 

Table  7-6:  Results  of  the  simulations  for  different  values  of  the  kinetic  constant,  diamond 

pillars,  variable  spacing. 

Settling  Time  for  k”: 

12  s 

Kinetic 

constant 

fc"  [m/s] 

Outlet  Propellant 

peroxide  mass  conversion 

fraction  [-]  [%] 

Temperature 

[K] 

Evaporation 
fraction  [%] 

le-3 

0.67 

23 

404 

72 

4.2e-3 

0.59 

33 

411 

91 

Table  7-7:  Results  of  the  simulations  for  different  values  of  the  kinetic  constant,  square  pillars, 

variable  spacing. 


In  general,  the  effects  of  changes  in  geometry  and  pillar  spacing  are  rather  negligible. 


7.6.6  Transient  response  of  the  device:  thermal-chemical  coupling 

The  successful  operation  of  the  device  relies  on  its  capability  of  producing  a  hot  gaseous  stream  at  its 
outlet. 


5  Reaction  rate  is  approximately  proportional  to  the  product  S/V  ■  CHz0z 
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This  in  turn  does  not  depends  solely  on  the  complete  decomposition  of  the  peroxide,  but  also  on  a  careful 
thermal  management  of  the  device  itself.  It  is  clear  that  there  is  competition  between  the  heat  generated 
inside  the  device  and  the  heat  losses  through  the  silicon  walls  towards  the  ambient.  The  ratio  of  heat  losses 
(a  surface  phenomenon)  to  heat  of  reaction  (a  volumetric  quantity)  increases,  as  the  device  is  downscaled, 
with  the  inverse  of  the  characteristic  dimension.  Furthermore,  due  to  the  high  thermal  conductivity  of  the 
silicon  substrate,  heat  transport  does  not  act  only  on  the  transversal  (towards  ambient)  direction,  but  has 
also  a  significant  component  on  the  stream-wise  direction.  This  means  that  a  plug-fow  1-D  analysis  may  fail 
considerably  in  predicting  the  outlet  flow  conditions,  since  the  usual  assumption  of  a  vaporization  process 
taking  place  in  a  thin  spatial  region  is  no  more  applicable.  Therefore,  also  the  fluid  vaporization  will 
somehow  tend  to  move  upstream. 

On  one  side,  this  can  be  thought  as  a  beneficial  effect  since  it  enhances  vaporization.  On  the  other  side, 
since  a  considerable  fraction  of  hydrogen  peroxide  vaporize  before  decomposing,  the  propellant 
concentration  decreases  (due  to  the  lower  mixture  density),  which  results,  according  to  the  assumed  1st 
order  reaction  kinetics,  to  a  substantially  reduced  rate  of  decomposition.  Furthermore,  a  drawback  may 
arise  due  to  the  fact  that  the  mixture  density  decreases  as  a  result  of  vaporization,  and  its  average  speed 
increases.  This  may  turn  the  reactor  towards  a  diffusion  limited  operation  mode  (too  high  u/D)  which 
basically  means  that  the  peroxide  flow  is  convected  away  from  the  reaction  pillars  too  fast,  before  having 
the  chance  to  decompose  completely. 

A  practical  view  of  the  above  phenomenon  is  seen  by  looking  at  the  temporal  evolution  of  the  outlet 
decomposition  in  the  example  below,  dealing  with  a  transient  simulation  performed  at  k"  =  4.5e-3  m/s, 
starting  from  a  converged  solution  obtained  at  k"  =  le-3  m/s  (i.e.  a  Mode#2  simulation).  Three  pictures  are 
shown,  depicting  the  evolution  of  the  average  outlet  a)  peroxide  concentration,  b)  vaporization  fraction  and 
c)  temperature  in  the  first  0.35  seconds  of  simulation  (temporal  resolution  is  0.05  s). 
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Figure  7-8:  Average  outlet  evolution  of  mass  fraction  (top),  evaporation  fraction  (middle)  and  temperature  (bottom). 

The  first  picture  is  representative  of  the  strong  coupling  between  thermal  effects  and  reaction  kinetics:  at 
the  first  time  step  the  outlet  peroxide  fraction  decreases  strongly  from  the  equilibrium  value  of  0.6  for  a  k 
=  le-3  m/s  (the  initial  condition)  down  to  0.44  due  to  the  enhanced  kinetic  k  =  4.5e-3  m/s.  At  subsequent 
times,  the  trend  is  reversed,  with  the  outlet  YH2q2  increasing  due  to  the  reduced  peroxide  concentration, 
which  is  in  turn  a  direct  consequence  of  the  ongoing  evaporation  (mixture  density  decreasing). 

On  one  side,  we  see  the  expected  behavior  for  both  the  temperature  and  vaporization  fraction,  which 
increase  monotonically.  On  the  other  side,  the  degree  of  conversion  at  the  outlet  is  only  marginally  higher 
than  the  one  obtained  at  the  lower  k  . 

A  similar  behavior  was  encountered  also  with  analysis  run  in  Mode#3,  for  sufficiently  short  raising  times  of 
the  kinetic  constant  (Table  7-6).  As  an  example,  in  Figure  7-9  results  are  shown  obtained  with  a  simulation 
where  k  was  allowed  to  raise  from  le-7  to  4.5e-3  m/s  in  6  seconds.  Again,  while  the  temperature  raises 
steadily,  the  outlet  concentration  does  not  decrease  in  a  monotonic  way,  with  a  bump  occurring  when  the 
fluid  is  evaporating. 

These  results  suggest  that  it  is  unlikely  that  a  complete  decomposition  will  be  reached  unless  much  higher 
temperatures  are  achieved  through  other  means  (externally  supply  through  e.g.  joule  heating).  This  later 
conclusion  however,  shall  be  further  investigated,  and  definitive  conclusions  can  be  drawn  once  a  full 
simulation  up  to  complete  evaporation  will  be  achieved. 

Finally,  it  is  worth  to  be  noted  that,  when  increasing  the  settling  time  for  the  kinetic  constant  up  to  12-14  s 
(as  per  Table  7-4,  Table  7-5  and  Table  7-7),  the  propellant  conversion  follows  instead  a  monotonic 
decrease. 
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Figure  7-9:  Prescribed  kinetic  constant  (top  panel)  vs  temperature,  evaporation  fraction,  and  outlet  peroxide 
concentration  obtained  with  a  Mode#3  simulation. 


7.6.7  Post  vaporization  phase 

Since  complete  vaporization  was  not  reached  with  a  single  simulation,  a  separated  treatment  of  the  post 
vaporization  phase  has  been  set  up:  this  consisted  in  using  the  same  geometry,  but  feeding  as  inlet  a 
gaseous  stream  with  a  temperature  equal  to  150°C  (i.e.  the  top  end  of  the  vaporization  interval)  and  a 
peroxide  mass  fraction  extrapolated  from  previous  results.  In  average,  the  highest  degree  of  vaporization 
was  achieved  for  a  peroxide  concentration  of  59%,  which  was  taken  as  input  for  the  present  simulations. 
With  this  single-phase  setup,  convergence  was  easily  achieved  with  steady  state  simulations. 

As  a  further  investigation,  we  checked  the  effect  of  the  change  in  order  of  magnitude  of  the  diffusivity 
coefficient,  whose  baseline  value  10'6  m2/s,  is  quite  high  for  the  gaseous  phase.  In  the  table  below,  the 
results  of  this  test  for  the  square  diamond  are  summarized,  for  values  of  the  kinetic  constant  k"  taken  as 
10'2  and  10 1  m/s. 
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Diffusivity  coeff. 

Kinetic  constant 

Outlet  peroxide 

Propellant 

Temperature 

[m2/ s] 

K  [m/s] 

mass  fraction  [-] 

conversion  [%] 

[K] 

le-2 

0.42 

52 

436 

le-8 

le-1 

0.15 

83 

463 

le-2 

0.42 

52 

436 

le-7 

le-1 

0.14 

84 

464 

le-2 

0.36 

54 

437 

le-6 

le-1 

0.10 

89 

464 

Table  7-8:  Results  for  the  post-vaporization  simulations 


The  degree  of  propellant  conversion  at  the  outlet  lies  in  between  83%  and  89%,  with  the  highest  value 
achieved  when  diffusion  coefficient  is  set  at  a  maximum,  as  expected. 


7.6.8  Summary 

The  results  of  the  post-vaporization  suggest  that  a  length  on  the  order  of  the  entire  modeled  domain  is 
needed  to  obtain  an  almost  complete  decomposition,  or  even  larger  as  the  diffusivity  reduces.  One  could 
then  be  tempted  to  assume  an  overall  length  being  twice  the  nominal  length:  one  to  achieve  full 
evaporation,  and  one  to  achieve  full  decomposition  within  the  gaseous  phase.  This  would  lead,  however,  to 
an  overestimation,  since  only  by  simulating  at  once  from  liquid  inlet  to  vapor  outlet  we  can  model 
consistently  the  thermal  coupling  of  the  substrate. 

7.7  A  step  forward  in  the  FEM  analysis 

We  know  that  the  decomposition  of  the  propellant  starts  as  soon  as  the  hydrogen  peroxide  gets  in  contact 
with  the  catalyst.  The  initial  temperature  is  usually  well  below  the  boiling  point  such  that  the  formed  water 
is  in  the  liquid  phase.  The  oxygen  is  in  the  gas  phase  and  it  appears  as  bubbles  in  the  liquid.  Since  the 
reaction  is  exothermic,  the  generated  heat  causes  a  rise  in  temperature  until  the  boiling  point  of  the  liquid 
mixture  is  reached.  The  increase  in  temperature  also  leads  to  an  increase  in  the  reaction  rate  leading  to 
more  gas  bubbles  in  the  liquid.  Some  bubbles  will  coalesce  with  other  bubbles  forming  larger  bubbles  or 
slugs.  Together  with  the  liquid  water,  part  of  the  peroxide  evaporates  before  it  can  decompose  in  the  liquid 
phase  leading  to  even  more  gas.  At  this  point,  rather  than  speaking  about  a  liquid  phase  with  gas  bubbles 
dispersed  in  it,  it  makes  more  sense  to  speak  about  a  gas  phase  with  liquid  droplets  dispersed  in  it. 
Eventually,  all  of  the  liquid  has  vanished  either  by  decomposition  or  evaporation.  Gaseous  peroxide  will 
decompose  as  resulting  in  a  further  increase  in  gas  temperature.  What  remains  is  a  gas  mixture  of  steam 
and  oxygen. 

Clearly,  all  the  features  above  cannot  be  captured  by  the  FEM  model  described  previously,  which  relies  on 
the  homogeneous  fluid  assumption:  an  advancement  in  the  level  of  detail  of  the  FEM  analysis  can  be 
pursued,  which  accounts  for  the  coexistence  within  the  dissociation  chamber  of  liquid  and  gas  as  separated 
phases.  The  first  steps  towards  such  development  is  outlined  in  the  following. 
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7.7.1  Modelling  Non-homogeneous  multiphase  flow  with  phase  change 

The  development  of  a  complete  computational  model  that  takes  into  account  the  entire  decomposition 
process  is  extremely  challenging.  Chemical  and  thermal  point  of  view  aside,  treatment  of  the  multiphase 
flow  component  during  the  initial  stages  of  the  decomposition  is  especially  difficult  and,  at  the  same  time, 
the  literature  is  very  limited. 

In  this  section,  we  report  the  groundwork  for  an  upcoming,  comprehensive  approach  to  the  multiphasic 
modelling  of  the  problem.  The  first  step  is  the  assessment  of  a  numerical  model  for  describing  phase 
change  in  the  proximity  of  a  solid  domain.  This  is  representative  of  the  catalytic  decomposition  at  the  pillar 
where,  other  than  the  chemical  reaction,  a  change  of  phase  occurs  (from  liquid  peroxide  to  gaseous 
oxygen). 

COMSOL  Multiphysics,  offers  several  physics  interface  to  simulate  different  kind  of  multiphase  flow 
problems,  therefore  initially  an  assessment  was  performed  to  select  the  most  suitable  for  the  present 
investigation.  Given  our  focus  and  in  particular,  the  need  to  capture  the  movements  of  the  fluid  interface 
and  of  incorporate  subsequently  COMSOL  physical  nodes,  as  the  chemical  reaction  and  thermal  ones,  we 
opted  for  the  Two-Phase  Level  Set  Interface.  This  physics  interface  solves  the  Navier-Stokes  equations 
(conservation  of  momentum)  and  a  continuity  equation  (conservation  of  mass),  according  to: 


p  —  =  V  ■  [-pi  +  p(Vu  +  Vur)]  +  Fg  +  Fst 
V  ■  u  =  0 


(7-19) 


Where  u  is  the  velocity  field,  p  the  pressure  field,  p  the  dynamic  viscosity  of  the  mixture,  Fg  the  gravity 
force  and  Fst  the  surface  tension  force  acting  on  the  interface  between  the  two  fluids. 

The  interface  position  is  tracked  by  solving  a  transport  equation  for  the  level  set  function: 


(7-20) 


The  left-hand  side  represent  interface  advection,  while  the  right-hand  side  term  includes  stabilization  and 
reinitialization  terms.  In  particular,  s  controls  the  interface  thickness  and  parameter  y  defines  the  intensity 
of  y  is  the  reinitialization  parameter  (set  to  1  by  default),  4>  is  the  level  set  function  (set  to  0.5  at  the 
interface)  and  £  is  the  interface  thickness  controlling  parameter  (set  as  default  to  hmax/2  where  hmax  is  the 
maximum  element  size  in  the  component).  The  level  set  function  is  fundamental  to  define  the  boundary 
between  the  immiscible  phases,  which  assumed  to  be  located  at  4>  =  0.5.  In  the  future,  the  level  set 
function  will  be  used  to  set  specific  regions  in  the  domain  where  interface  phenomena  occur,  such  as 
evaporation. 

The  main  objective  is  to  model  the  production  of  the  gas  phase  at  the  catalytic  surface.  For  this  purpose,  a 
solid  theoretical  base  arrives  from  studies  related  to  phase  change  phenomena  [7]  .  They  are  very  useful,  in 
our  case,  for  two  reasons: 

•  We  can  approximate  the  production  of  oxygen  as  a  phase  change  in  a  small  region  around  the 
surface  of  the  pillars; 

•  We  can  model  the  evaporation  of  water  and  peroxide  due  the  exothermic  chemical  reaction. 
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The  governing  equations  involved  to  simulate  an  isothermal,  laminar  and  incompressible  two-phase  flow 
with  a  phase  change  are  essential  the  default  relations,  seen  previously,  with  the  addition  of  some  source 
terms  both  in  the  continuity  and  in  the  level  set  transport  equations: 


V  ■  u  =  m 


^  +  V  •  («4>)  =  yV  ■  (sV4>  -  4>(1  -  'W  j^jj  Pl 


1  1 

-Pv  Pl- 


Vcj)  \  m 


(7-21) 

(7-22) 


where  m  is  the  rate  of  phase  change,  pt  and  pv  are  the  density  of  the  liquid  and  gas  phase,  respectively. 
The  source  terms  have  been  implemented  in  COMSOL  Multiphysics  as  weak  contributions. 

This  method  models  the  change  of  phase  as  a  volume  phenomenon,  while  ideally,  we  would  like  it  to 
happen  as  a  boundary  one,  i.e.  occurring  at  the  reactive  surface.  This  drawback  was  overcome  by  confining 
the  source  terms  within  a  small  volume  surrounding  the  pillar  (see  Figure  7-11). 

Finally  we  underline  also  that  the  use  of  the  conservative  form  in  Eq.  (7-22)  was  found  to  be  of  paramount 
importance  to  ensure  mass  conservation:  by  using  conservative  form,  the  level  set  function  assumes  the 
meaning  of  volume  fraction,  rather  than  that  of  an  advected  physical  boundary,  as  in  Eq.  (7-20). 

7.7.2  A  mass  conservation  test 

To  validate  the  implemented  model  a  simple  conservation  test  was  set  up.  The  purpose  of  which  is  to  verify 
that: 

•  the  amount  of  gas  generated  is  consistent  with  the  prescribed  rate  of  production  m; 

•  the  mass  of  gas  generated  corresponds  to  the  mass  of  liquid  consumed  due  to  phase  change. 

To  this  end,  we  considered  a  simplified  domain  with  a  single  diamond  pillar  (Figure  7-10).  Apart  from  the 
outlet  on  the  right,  all  the  boundaries  are  considered  as  wall  with  no-slip  condition.  As  initial  value,  the  blue 
region  represents  the  liquid  domain  while  the  red  region  represents  the  gas  domain;  the  latter  is  included 
to  avoid  that  liquid  phase  getting  out  of  the  domain,  therefore  ensuring  that  liquid  phase  disappears  only 
due  to  phase  change. 


Time=0  s  Surface:  Volume  fraction  of  fluid  1  (1) 
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Figure  7-10:  Volume  fraction  of  gas  at  time  0s 


Figure  7-11:  Geometric  domain  with  highlighted  the  region 
where  gas  phase  creation  occurs 
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Results  are  shown  in  the  following,  in  particular: 

•  Figure  7-12  that  the  mass  of  gas  produced  during  time  (blue  line)  reaches  the  theoretical  value 
prescribed  by  the  chosen  source  m  (green  line); 

•  Figure  7-13  shows  that,  parallel  to  gas  production,  an  equal  amount  of  liquid  phase  exits  the 
domain. 

As  stated  also  at  the  end  of  the  previous  section,  the  mass  conservation  is  not  satisfied  when  implementing 
interface  advection  in  non-conservative  form. 


Figure  7-12:  Mass  of  gas  included  in  the  domain  Figure  7-13:  Mass  of  liquid  included  in  the  domain 
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Figure  7-14:  Snapshots  of  the  temporal  evolution  of  gaseous  phase  created  at  the  pillar  during  the  mass  conservation 

test. 


Page  33  of  37 


DISTRIBUTION  A.  Approved  for  public  release:  distribution  unlimited 


Final  Report 


Version 

1.0 

Date 

28/02/2017 

8.  MICROTHRUSTER  DESIGN 

The  numerical  analysis  performed  so  far  provided  some  good  insight  into  the  physical  mechanisms 
governing  the  decomposition  of  the  propellant.  Nevertheless,  quantitative  guidelines  for  an  optimal  design 
of  the  chamber  cannot  be  drawn  with  enough  confidence  such  to  select  a  priori  a  single  preferred  design 
for  prototyping. 

Furthermore,  manufacturing  constraints  needs  to  be  considered:  for  example,  despite  one  would  like 
having  the  catalytic  surfaces  as  close  as  possible  (i.e.  low  pillar  spacing),  on  the  other  side,  this  would  result 
in  decreased  accessibility  of  the  pillars  lateral  surface  when  it  comes  to  catalyst  deposition.  Catalyst, 
whether  platinum  or  palladium,  is  deposited  through  PVD  process  using  an  e-gun  in  a  vacuum  chamber 
equipped  with  a  planetary  sample  holder  undergoing  a  "wobbling"  rotational  motion.  The  uniformity  of 
the  deposition  is  difficult  to  be  addressed  quantitatively:  this  can  be  done,  partially,  only  a-posteriori,  by 
cutting  the  devices  and  inspecting  the  resulting  sections  with  a  scanning  electron  microscope. 

For  this  reasons,  for  the  design  of  the  next  generation  of  prototypes  an  approach  was  chosen  which  tries  to 
explore  a  design  space  as  much  as  possible:  the  main  design  parameters  are  the  length,  width,  shape/size 
of  the  pillars  and  spacing.  This  approach  is  made  possible  since  the  manufacturing  process  (DRIE,  deep 
reactive  ion  etching),  allows  to  accommodate  tens  of  devices  within  a  single  process:  etching  on  a  6" 
diameter,  625  pm  thick  silicon  wafers. 

Two  pillars  geometry  were  considered:  diamond  and  arrow-like,  and  for  each  of  them  3  different  widths,  2 
lengths  and  whether  uniform  or  variable  spacing  pillars  have  been  designed.  This  results  in  a  total  of  24 
geometries  which  are  due  for  production.  A  sample  drawing  for  one  of  this  geometry  is  shown  in  Figure 
8-1. 


Figure  8-1:  Example  of  1  out  of  24  thruster  geometries  due  for  manufacturing. 


8.1  Design  of  a  test  bench  for  the  microthruster 

Testing  of  microfluidic  devices  is  quite  a  complex  task,  and  often  the  design  of  an  ad-hoc  test  bench  is 
mandatory  for  an  efficient  handing  of  the  experimental  activities. 

The  microthruster  is  essentially  composed  by  an  inlet  region  with  a  circular  hole  at  the  base,  a 
decomposition  chamber  and  an  exit  region  with  a  convergent/divergent  nozzle.  At  the  base  of  the  central 
region,  characterized  by  a  series  of  catalytic  pillars,  there  are  applied  a  series  of  electrical  resistances  with 
two  purposes:  heating  this  region  to  enhance  the  rate  of  reaction;  measuring  the  temperature  as  resistance 
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thermometer.  Accordingly,  the  test  bench  (Figure  8-2)  must  enable  both  a  fluid  dynamic  connection,  for  the 
propellant  feeding,  and  an  electrical  connection  for  the  wiring  system. 


Figure  8-2:  Concept  for  the  test  bench 

For  the  preliminary  design  of  a  test  bench  there  are  some  problems: 

•  Centering  of  the  device; 

•  Leakage  at  the  inlet; 

•  Coupling  between  the  fluid  dynamic  connection  and  the  electrical  one. 

The  concept,  analyzed  in  this  section,  tries  to  solve  these  issues  with  a  solution  that  is  reliable,  practical  and 
flexible  for  different  size  of  microthrusters. 

About  the  first  point,  considering  the  flexibility  feature  of  the  system,  it  is  necessary  to  center  the  device 
through  the  inlet  by  a  cylindrical  component  (Figure  8-3). 


Figure  8-3  Lower  part  of  the  fluid  dynamic  connection  Figure  8-4:  Exploded  diagram  of  the  fluid  dynamic  connection 
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The  seal  at  the  entry  hole  of  the  propellant  is  guaranteed  by  two  gaskets  (Figure  8-4),  pressed  against  the 
thruster  through  a  bolt-nut  system. 

The  electrical  connection,  is  essentially  realized  through  a  series  of  spring  contact  probes.  For  this  to  work  it 
is  fundamental  using  a  non-conductive  material  and  an  accurate  centering  to  ensure  an  effective  contact 
between  the  probes  and  the  heating/sensing  elements. 


9.  CONCLUSIONS 

This  report  discusses  the  activities  performed  under  USAF  AFOSR  Grant  n°  FA9550-16-1-0081.  These  have 
been  part  of  a  wider  collaboration  of  the  team  within  a  consortium  for  the  development  of  HTP  MEMS- 
based  micro-thrusters  in  the  mN  range  thrust. 

The  most  relevant  contribution  is  the  development  of  a  numerical  model  which  jointly  simulates  the 
chemical,  fluid-dynamics  and  thermal  behavior  of  the  propellant  stream  inside  the  combustion  chamber,  as 
well  as  the  thermal  coupling  with  the  silicon  substrate.  Within  this  model,  an  element  of  novelty  is 
represented  by  the  inclusion  of  the  evaporation  equations  based  on  the  variable  heat  capacity  method  and 
the  concept  of  pseudo-phase.  This  method,  offers  a  tractable,  yet  consistent,  way  to  model  the  evaporation 
in  a  chemical  reacting  mixture.  As  a  main  result  of  the  analysis,  it  emerged  that  the  silicon  substrate  has  a 
great  impact  on  the  flow  evolution,  such  that  the  evaporation  region  is  distributed  all  over  the  length  of  the 
device.  Transient  simulations  allowed  to  get  more  insight  into  the  coupling  of  thermal  and  chemical  related 
phenomena,  showing  a  non-monotone  degree  of  conversion  at  the  outlet  due  to  the  decrease  of  the 
propellant  density. 

A  sensitivity  study  was  performed  to  insect  the  impact  of  variable  geometry  and  non-equispaced  pillars. 
Overall,  the  impact  on  the  degree  of  conversion  achieved  at  the  outlet  was  found  to  be  quite  low. 

A  strategy  was  also  outlined  for  increasing  further  the  level  of  detail  of  the  simulation  with  a  non- 
homogeneous  multi-phase  analysis,  based  on  a  modification  of  the  level-set  method  to  allow  for  the 
change  of  phase  at  the  catalytic  surface.  The  consistency  of  the  approach  has  been  validated  through  mass 
conservation  tests,  with  good  results.  Addition  of  chemical  reactions  and  thermal  equations  is  till  to  be 
performed.  The  high  computational  burden  suggests  that  such  method  hardly  could  be  extended  to  include 
a  fully  coupled  multiphysics  simulation  on  a  3D  domain:  most  likely,  this  approach  will  rather  be  confined  to 
more  fundamental  research  objectives,  such  as  validation  of  the  homogeneous  model  in  a  simplified,  lower 
dimensional  test  case. 

Some  design  activities  were  also  performed  in  preparation  of  the  manufacturing  of  the  next  batch  of 
prototypes.  With  respect  to  the  previous  ones,  these  will  incorporate  actuators  and  sensors  in  the  form  of 
heating  elements  and  resistive  temperature  sensors.  Different  combustion  chambers  were  designed,  with 
equispaced  or  variable  spaced  pillars. 

Since  the  past  experience  with  the  previous  prototypes  showed  that  the  interfaces  are  a  critical  part  of  the 
system  design,  a  dedicated  mounting  support  has  been  designed  which  will  allow  both  fluidic  and  electrical 
reliable  connections. 
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The  batch  of  prototypes  are  due  to  production  at  an  external  research  institution  in  the  next  weeks.  Results 
of  both  the  experimental  activities  and  of  the  numerical  ones  will  be  published  on  the  open  literature,  with 
the  due  acknowledgement  of  the  funding  source. 
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